######################################################################
# Time-stamp: <05/08/20 09:44:40 ostolop>
#
# Package: EP::EPC::Core::Threeway_Similarity_Analysis
#
# Three-way Similarity Analysis EPC
######################################################################
package EP::EPC::Core::Threeway_Similarity_Analysis;

use strict;

use EP::Config;

use EP::Common;
use EP::Common::General;
use EP::Common::Error qw ( :try );

use EP::EPC::Core;
use EP::EPC::Core::Clustering::exp;

use EP::Visualization::Grid;
use EP::Visualization::ThreewaySimilarityGraphic;

use Bio::DB::GFF;
use Bio::Graphics;
use Bio::SeqFeature::Generic;

use File::Temp;

use DBI;

use PDL;
use PDL::NiceSlice;

# use SVG;

use GO::TermFinder;
use GO::AnnotationProvider::AnnotationParser;
use GO::OntologyProvider::OntologyParser;

use base qw ( EP::EPC::Core::Component );


# DB Connection Info 
use constant PG_DBI_DRIVER => "dbi:Pg:dbname=";
use constant PG_ADAPTOR => "dbi::pg";
use constant DBNAME => "gff";
use constant USER => "postgres";
use constant PWD => "";

use constant SEQUENCE_ALIGN_DIR => "$EP::Config::h->{TSA_DATA_PATH}/align/";
use constant INTERPRO_MATCH_DIR => "$EP::Config::h->{TSA_DATA_PATH}/interpro/result/";


sub process_component {

  my $self = shift;

  my ($svg_file, $width, $height) = $self->run_threeway_similarity_analysis ();

  my $process_result = {
			visualization  => "$EP::Config::h->{HTML_HREF}/$EP::Config::h->{WWW_TMP}/$svg_file",
			visualization_width => $width,
			visualization_height => $height,
			src_dataset_id => $self->{src_dataset}->id,
			dst_dataset_id => $self->{src_dataset}->id
		       };

  $self->{CD}->{process_result} = $process_result;

  return OK;
}

sub run_threeway_similarity_analysis {

  my $self = shift;

  my $dataset = $self->{src_dataset};
  my $input_gene =  $self->{query}->param ( "input_gene" );
  my $input_gene_id = EP::EPC::Core::Clustering::exp::get_gene_ids ($dataset, $input_gene); 
  my $input_gene_name = $dataset->{data}->{row}->{$input_gene_id->[0]}->{full_annot}->{NAMES};
 
  chomp($input_gene);

  # computing the similarity  
  $self->{session}->status ( "Computing Expression Similarity..." );
  my $expression_similarity = $self->_expression_similarity ($input_gene, $input_gene_id, $dataset);

  $self->{session}->status ( "Computing Sequence Similarity..." );
  my $sequence_similarity = $self->_sequence_similarity ($input_gene_id->[0], $input_gene_name, $dataset);

  $self->{session}->status ( "Computing Motif Similarity..." );
  my $motif_similarity = $self->_motif_similarity ($input_gene_id->[0], $input_gene_name, $dataset);


  # initializing Gene Ontology
  $self->{session}->status ( "Initializing Gene Ontology..." );
  $self->_go_initialize();


  # creating the data for the visualization of the result 
  $self->{session}->status ( "Creating the data for visualization..." );
  my $data = $self->_create_data (
                     $expression_similarity, 
                     $sequence_similarity, 
                     $motif_similarity,
                     $input_gene_name );
 
  # debugging
  # &_dump( "/usr/local/EPNG/tmp/debug.txt", $data, "data_debugging" );

  # creating the visualization file 
  $self->{session}->status ( "Creating the visualization file..." );
  my ($svg_file, $img_width, $img_height) = $self->create_visualization ($data);

  
  return  ($svg_file, $img_width, $img_height);

}

# Sub: _expression_similarity
#
# PURPOSE:  Compute the expression similarity between an input gene and all other pair genes 
# ARG:      $input_gene      ex) YAL038W
#           $input_gene_id   ex) index no of input gene in the expression dataset  
#           $dataset         ex) dataset contains the expression dataset
# RETURN:   $expression_similarity
sub _expression_similarity {

  my $self = shift;
  my $input_gene = shift;
  my $input_gene_id = shift;
  my $dataset = shift;

  my $distance_measure = $self->{query}->param ( "distance" );
  my $gene_no = $dataset->data->raw_data->getdim(1);

  EP::EPC::Core::Clustering::exp::create_similarity_search( 
                 $dataset, 
                 $distance_measure, 
                 $input_gene_id, 
                 $gene_no, 
                 1 ); # 1 to impute missing values ???  => This could make the inaccurate result!

  # $dataset->{clustering} contains the filename of the distance matrix.
  die "Couldn't calculate distances!" if not $dataset->{clustering};

  my $clustering_filename = "$EP::Config::h->{EPNG_USRDATA_PATH}/" . 
			     $dataset->folder->location . 
                            "/$dataset->{clustering}"; 

  # parse the result clustering file and return the similarity result 
  return $self->_parse_clustering_file($clustering_filename);

}

# Sub: _parse_clustering_file
#   
# PURPOSE: parse the result clustering file from CLUSTERING module and convert it into expression similarity 
sub _parse_clustering_file {
 
  my $self = shift;
  my $filepath = shift;
  my ( @members, @distances );

  open CLUST, $filepath or die "Couldn't open clustering results file $filepath !";

  while (<CLUST>) {
    chomp;
    my ($cluster_nr, $cluster_size) = $_ =~ /CLUSTER (\d+?) SIZE=(\d+?)/ if /CLUSTER \d+/;
    @members = split /\s+/ if /MEMBERS:/;
    @distances = split /\s+/ if /DISTANCES:/;
  }

  close CLUST;

  # shift - because the first item is "MEMBERS:" or "DISTANCES:"
  shift @members;
  shift @distances;

  # (jaak's code returns 1-based ID's)
  $_-- foreach @members;

  map { chomp } @members;
  map { chomp } @distances;

  # key should be the gname not the expression source dataset id
  for(my $i=0; $i<=$#members; $i++) {
    my $gname = $self->{src_dataset}->{data}->{row}->{$members[$i]}->{full_annot}->{NAMES};
    $members[$i] = $gname;
  }

  # convert distance to similarity
  $self->_distance2similarity( \@distances, $distances[$#distances], $#distances+1 ); 

  # currently input gene is selected from the expression dataset,
  # therefore, its data is certainly valid
  my %result;
  $result{valid} = 1;
  @result{@members} = @distances;
  return \%result;

}

# Sub: _distance2similarity
#   
# PURPOSE: convert distance measure into similarity measure in the given data
sub _distance2similarity {

  my $self = shift;
  my $distance = shift;
  my $max_distance = shift; 
  my $gene_no = shift;

  # set the similarity of the most dissimilar gene to 0
  # set the similarity of the input gene itself to 1
  # other genes has the relative value 
  for(my $i=0; $i<$gene_no; $i++) {
	$distance->[$i] = ( $max_distance - $distance->[$i] ) / $max_distance;
  }  

}

# Sub: _sequence_similarity
#
# PURPOSE:  Compute the sequence similarity between an input gene and all other pair genes 
# ARG:      $input_gene      ex) YAL038W
#           $input_gene_id   ex) index no of input gene in the expression dataset  
#           $dataset         ex) dataset contains the expression dataset
# RETURN:   $sequence_similarity
sub _sequence_similarity {

  my $self = shift;
  my $input_gene_id = shift;
  my $input_gene_name = shift; 
  my $dataset = shift; 

  my %result;
  my $input_gene_seq_filename = SEQUENCE_ALIGN_DIR . $input_gene_name . ".summary";

  # find if $input_gene has the sequence similarity information
  # if not, return with the tag this info is invalid !  
  $result{valid} = 0;
  return \%result if not -e $input_gene_seq_filename; 
  	
  # if it has, return with the sequence similarity info
  $result{valid} = 1;
  return &_parse_sequence_file($input_gene_seq_filename, \%result);

}

# Sub: _parse_sequence_file
#
# PURPOSE:  parse the precomputed sequence similarity file between an input gene and all other pair genes 
sub _parse_sequence_file {

  my $filepath = shift;
  my $result = shift;
  my @line;

  open FH, $filepath or die "Couldn't open the sequence align result file $filepath !";

  while (<FH>) {
    chomp;
    @line = split /\t+/;

    # $line[0] = ORF, $line[1] = sequence similarity
    $result->{$line[0]} = $line[1];
  }

  close CLUST;

  return $result;
}


# Sub: _motif_similarity
#
# PURPOSE:  Compute the motif similarity between an input gene and all other pair genes 
# ARG:      $input_gene_id   ex) index no of input gene in the expression dataset 
#           $input_gene_name ex) YAL038W
#           $dataset         ex) dataset contains the expression dataset
# RETURN:   $expression_similarity
sub _motif_similarity {

  my $self = shift;
  my $input_gene_id = shift;
  my $input_gene_name = shift;
  my $dataset = shift; 

  my %result;

  my $input_gene_motif_filename = INTERPRO_MATCH_DIR . $input_gene_name . ".summary";

  # find if $input_gene has the sequence similarity information
  # if not, return with the tag this info is invalid !  
  $result{valid} = 0;
  return \%result if not -e $input_gene_motif_filename; 
  	
  # if it has, return with the sequence similarity info
  $result{valid} = 1;
  return &_parse_motif_file($input_gene_motif_filename, \%result);

}

# Sub: _parse_motif_file
#   
# PURPOSE:  parse the precomputed motif similarity file between an input gene and all other pair genes 
sub _parse_motif_file {

  my $filepath = shift;
  my $result = shift;
  my @line;

  open FH, $filepath or die "Couldn't open the interpro motif result file $filepath !";

  while (<FH>) {
    chomp;
    @line = split /\t+/;

    # $line[0] = ORF, $line[1] = interpro match 
    $result->{$line[0]} = $line[1];
  }

  close CLUST;

  return $result;
}


# Sub: _go_initialize
#
# PURPOSE:  Initialize the GO Ontology Provider for extracting GO terms for each gene 
sub _go_initialize {

    my $self = shift;
    my $species = "scerevisiae";

    # obtain the process ontology structure file
    $self->{session}->status ( "Reading Process Ontology Tree ..." );
    my $process_ontology_file = "$EP::Config::h->{EPNG_GODATA_PATH}/process.ontology";
    if ( -f $process_ontology_file ) {
	$self->{ontologyProvider}->{P} = GO::OntologyProvider::OntologyParser->new( ontologyFile => $process_ontology_file );
    } else {
	die "Could not find ontology file for $process_ontology_file $!";
    }

    # obtain the function ontology structure file
    my $function_ontology_file = "$EP::Config::h->{EPNG_GODATA_PATH}/function.ontology";
    if ( -f $function_ontology_file ) {
	$self->{ontologyProvider}->{F} = GO::OntologyProvider::OntologyParser->new( ontologyFile => $function_ontology_file );
    } else {
	die "Could not find ontology file for $function_ontology_file $!";
    }
    
    # obtain the mapping file for this species
    $self->{session}->status ( "Reading Annotation Mappings ..." );
    my $annotation_file = $EP::Config::h->{EPNG_GODATA_PATH} . "/" . $EP::Config::h->{SPECIES_METADATA}->{$species}->{annotationFile};
    if ( -f $annotation_file ) {
	$self->{annotationProvider} = GO::AnnotationProvider::AnnotationParser->new( annotationFile => $annotation_file );
    } else {
	die "Could not find the annotation file for species: $annotation_file";
    }
   

    # finding out the representative GO terms for this group of genes can be added here ! 
    #$self->{session}->status ( "Initializing GO TermFinder ..." );
    #$self->{termFinder}->{P} = GO::TermFinder->new( annotationProvider => $self->{annotationProvider},
    #					       ontologyProvider => $self->{ontologyProvider}->{P},
    #					       aspect => "P" );
    #$self->{termFinder}->{F} = GO::TermFinder->new( annotationProvider => $self->{annotationProvider},
    #					       ontologyProvider => $self->{ontologyProvider}->{F},
    #					       aspect => "F" );
}

# Sub _create_data
#
# PURPOSE:  create the data which contains all similarity information and go annotation  
# ARG:      $e_similarity   : expression similarity dataset
#           $s_similarity   : sequence similarity dataset
#           $m_similarity   : motif similarity dataset
#           $input_gene_name:
# RETURN:   $data            
sub _create_data { 
 
  my $self = shift; 
  my $e_similarity = shift;
  my $s_similarity = shift;
  my $m_similarity = shift;
  my $input_gene_name = shift; 

  my $data;
  my $gene_no=0;
  my $valid_gene_no=0;
 
  # init the genename and gid from GFF db and return the total number of genes
  my ($dbh, $sql, $sth, $result_no, $row, $gid, $gname);

  $dbh = DBI->connect ( $EP::Config::h->{SPECIES_METADATA}{yeast}{GFF_dsn}, "postgres", "", { RaiseError=>1, AutoCommit=>0 } );

  # currently only 'Y~~' ORF are selected
  $sql = "select gid, gname from fgroup where gname like \'Y%\' order by gid";
  $sth = $dbh->prepare($sql);
  $sth->execute();
                                   
  # previously to indicate the gene which doesn't have info, -1 is used.
  # but now we use zero as their similarity.  
  while($row = $sth->fetchrow_arrayref) {
        $gid = ${$row}[0];
        $gname = ${$row}[1];

        # first assume that this gene has 3 similarity measures
	$data->{$gname}->{valid} = 1;

        # EXPRESSION
        if( $e_similarity->{valid}==0 ) {
 		$data->{$gname}->{expression} = 0;   
		$data->{$gname}->{valid} = 0;
                $data->{$gname}->{exp_id} = -1; 
        }
        elsif( exists( $e_similarity->{$gname} ) ) {
		$data->{$gname}->{expression} = $e_similarity->{$gname};
        }
        else {
 		$data->{$gname}->{expression} = 0;   
		$data->{$gname}->{valid} = 0;
                $data->{$gname}->{exp_id} = -1; 
	}	

        # SEQUENCE
        if( $s_similarity->{valid}==0 ) {
 		$data->{$gname}->{sequence} = 0;   
		$data->{$gname}->{valid} = 0;
        }
        elsif( exists( $s_similarity->{$gname} ) ) {
		$data->{$gname}->{sequence} = $s_similarity->{$gname};
        }
        else {
 		$data->{$gname}->{sequence} = 0;     
		$data->{$gname}->{valid} = 0;
	}	

        # MOTIF 
        if( $m_similarity->{valid}==0 ) {
 		$data->{$gname}->{motif} = 0;   
		$data->{$gname}->{valid} = 0;
        }
        elsif( exists( $m_similarity->{$gname} ) ) {
		$data->{$gname}->{motif} = $m_similarity->{$gname};
        }
        else {
 		$data->{$gname}->{motif} = 0;        
		$data->{$gname}->{valid} = 0;
	}	


        if( $data->{$gname}->{valid} == 1 ) {
                $valid_gene_no++;

                # this part should go to the outside of this if statement
                # if we deal with all genes regardless of their validity which means that it haa all three data
                $self->_get_go_annotation($data, "P", $gname);
                $self->_get_go_annotation($data, "F", $gname);

                # store the exp_id for valid genes
                $data->{$gname}->{exp_id} = EP::EPC::Core::Clustering::exp::get_gene_ids ($self->{src_dataset}, $gname)->[0]; 
	}

        $gene_no++;

  }


  # remember how many gene name is included in $data
  $data->{META_valid_gene_no} = $valid_gene_no;
  $data->{META_gene_no} = $gene_no;
  $data->{META_exp_dataset} = $self->{src_dataset}->{data}->{matrix};
  $data->{META_input_gene} = $input_gene_name; 


  # load the motif image url from the precomputed file
  $self->_load_motif_url( $data );


  return $data;  
}

# Sub _get_go_annotation
#
# PURPOSE: get go annotation for each gene in the data 
sub _get_go_annotation {

  my $self = shift;
  my $data = shift;
  my $aspect = shift;
  my $gname = shift;

  # check if $gname is ambiguous
  if( $self->{annotationProvider}->nameIsAmbiguous($gname) ) {
      $data->{$gname}->{"go_$aspect"}->{AMBIGUOUS} = 1;
      return;
  }

  # first find the go ids by $gname
  my $goIdsRef = $self->{annotationProvider}->goIdsByName( name => $gname, aspect =>$aspect ); 

  # for each go id, get the ontology terms
  foreach my $goId ( @$goIdsRef ) {
	my $node = $self->{ontologyProvider}->{$aspect}->nodeFromId($goId);
     
        # register the go term for gname
        if( defined($node) ) {
            $data->{$gname}->{"go_$aspect"}->{$goId} = $node->term;
        }
        else {
            $data->{$gname}->{"go_$aspect"}->{$goId} = "not defined";
	    warn "[WARNING] no go term node is defined for go id ($aspect) $goId ";
        }
  }

}

# Sub _load_motif_url
#
# PURPOSE: load motif image url from the precomputed file 
sub _load_motif_url {

  my $self = shift;
  my $data = shift;

  my $filepath = "$EP::Config::h->{TSA_DATA_PATH}" . "/interpro/Protein_Motif_URL.txt"; 
  open FH, $filepath or die "Couldn't open clustering results file $filepath !";

  while (<FH>) {
    chomp;
    my @members = split /\t+/;
    if( exists($data->{$members[0]}) ) {
        $data->{$members[0]}->{motif_url} = $members[1];
        $data->{$members[0]}->{motif_no} = $members[2];
        $data->{$members[0]}->{motif_protein} = $members[3];
    }
  }

  close FH;
}



# Sub: _create_visualization
#
# PURPOSE:  Create the visualization file by calling the ThreewaySimilarityGraphic.pm
# ARG:      $data : generate from _create_data function
# RETURN:   $svg_filename
#           $svg_width
#           $svg_height
sub create_visualization {
  my $self = shift;
  my $data = shift;			

  my $grid = new EP::Visualization::Grid ( 1, 1 );
  my $tsa  = new EP::Visualization::ThreewaySimilarityGraphic ( data => $data );

  $grid->fillCellArray ( [ [ $tsa ] ] );

  # this $svg_file is generated for displaying purpose
  # the content of $svg_file comes from $grid
  my $svg_file = new File::Temp ( DIR => $EP::Config::h->{TMP},
				  UNLINK => 0,
				  SUFFIX => ".svg" );

  print $svg_file $grid->to_svg( zoom => 0 );
  close $svg_file;

  # only return the file name 
  # don't remove the ()
  my ($svg_filename) = $svg_file =~ /$EP::Config::h->{TMP}\/(.*)/;

  return ( $svg_filename,
	   $grid->getGridWidth(),
	   $grid->getGridHeight() );
}


1;
